Range and Velocity Estimation in Radar using Maximum Likelihood Estimator¶

Based on¶

  • Range Estimation in Radar using Maximum Likelihood Estimator, H. Sadia et.al.
  • Efficient Single Frequency Estimators, Vaughan Clarkson

Generating time-base for calculations¶

In [3]:
num_samples = 1500
sample_freq = 1000
sample_period = 1/sample_freq
# t is shifted so that the center is around 0 (ie goes from -x to +x) for pulse shaping
t = (np.arange(num_samples) - num_samples/2) * sample_period

# unshifted version:
t_noshift = np.arange(num_samples) * sample_period

Generating TX Pulse¶

  • Range Estimation in Radar using Maximum Likelihood Estimator uses a rectangular pulse for analysis, this is unrealistic (infinite BW)
  • We use a gaussian pulse, which is commonly used in radar eg.
    • A Gaussian pulse generator for Ultra-WideBand radar system, M. Dhieb et.al. 2007
    • Tunable Subnanosecond Gaussian Pulse Radar Transmitter: Theory and Analysis, R. Feghhi et.al. 2020
In [4]:
carrier_freq = 20

# How much will the pulse shaping reduce the bandwidth?:
pulse_fractional_bw = 0.2
x = gausspulse(t, carrier_freq, bw=pulse_fractional_bw)
In [5]:
fig = px.line(x=t, y=x, title="Transmitted Signal")
fig.show()

Generating a simulated RX pulse¶

  • Delay the timebase an arbitrary amount of seconds to simulate return delay
  • Shift the frequency an arbitrary amount to simulate doppler shift
  • Add some AWGN
In [43]:
rx_delay_time = 0.123
doppler_multiple = 1.123
rx_delay_samples = rx_delay_time * sample_freq
y = gausspulse(t - rx_delay_time, carrier_freq * doppler_multiple, bw=pulse_fractional_bw)
In [45]:
awgn = np.random.normal(0, 0.3,  num_samples)
y_awgn = y + awgn

Add AWGN to recieved signal (Fig. 5 in paper). It can tolerate large amounts of AWGN because AWGN has no autocorrelation, and thus will not impact the correlation result at all.

In [67]:
dataframe = [dict(t=t_iter, tx=tx_iter, rx=rx_iter) for t_iter, tx_iter, rx_iter in zip(t, x, y_awgn)]
fig = px.line(dataframe, x="t", y=["tx", "rx"], facet_row="variable", color=px.NO_COLOR, labels=dict(value="Amplitude", t="Time (s)"))
fig.show()

MLE Estimate of delay¶

TODO derive mle -> correlation here from paper¶

In [47]:
corr = np.correlate(x, y_awgn, mode="full")[:num_samples]
corr_t = np.arange(len(corr))

Correlation corr is equivalent to the $z[n]$ in the paper.

$z[n] = \sum_{n=n_0}^{n_0+M-1} y[n]\cos(\omega_0(n-n_0))$

This is equivalent to cross-correlation, as we are multiplying a signal $y[n]$ with another shifting signal $\cos(\omega_0(n-n_0))$

This should be equivalent to Fig. 6 and Fig. 7 in paper.

Unlike the paper, we are only taking the first half of the correlation $z[n]$. These parts correspond to positive delays. Since there is not situation in which the recieved pulse would come in BEFORE the tx pulse was transmitted (ie. a negative delay), these can be safely ignored.

In [48]:
max_corr = np.argmax(corr)
print("argmax z[n]:", max_corr)
print("Actual delay:", rx_delay_time)
print("Estimated delay:", t_noshift[num_samples - 1 - max_corr])
argmax z[n]: 1376
Actual delay: 0.123
Estimated delay: 0.123

In the paper, they recieve a result of max_corr = 100 - delay.

We get max_corr = 99 - delay. I am assuming this is due to MATLAB (used in paper) having 1-indexing??

In [49]:
fig = px.line(x=corr_t, y=corr, title="Correlation, z[n]")
fig.add_vline(x=max_corr, line_color="red")
fig.show()

MLE Estimate of doppler shift¶

TODO derive MLE -> PSD max here¶

In [68]:
psd_f, psd_x = periodogram(y_awgn, fs=sample_freq, nfft=4*len(y_awgn))
fig = px.line(x=psd_f, y=psd_x, title="PSD of rx signal")
fig.add_vline(x=psd_f[np.argmax(psd_x)], line_color="red")
fig.show()
In [51]:
est_rx_freq = psd_f[np.argmax(psd_x)]
est_doppler_multiple = est_rx_freq / carrier_freq

print("Actual doppler shift:", doppler_multiple)
print("Estimated doppler shift:", est_doppler_multiple)
Actual doppler shift: 1.123
Estimated doppler shift: 1.125

Based on Efficient Single Frequency Estimators, Vaughan Clarkson. Claims that the MLE of frequency of recieved signal is equal to the maximum of the PSD.

$\hat{\omega}_{MLE} = argmax_\theta\bigl\lvert \sum_{n=0}^{N-1} x_n e^{-j\theta n} \bigr\rvert$

Link to paper

In [ ]: